Descriptive Maps

library(rgdal)
library(ggplot2)
library(readr)
library(dplyr)

Function

plot_fort_maps2 <- function(data, fill_str, legend_title, plot_title, size=0.1){
  require(ggplot2)
  require(ggsn)
  require(broom)
  ggplot() +                                               # initialize ggplot object
    geom_polygon(                                          # make a polygon
      data = data,                                    # data frame
      aes_string(x = "long", y = "lat", group = "group",                # coordinates, and group them by polygons
                 fill = fill_str),
      size=size, color="black") +                # variable to use for filling
    scale_fill_brewer(name=legend_title, palette = "RdYlBu", direction = -1,
                      drop = FALSE) + # fill with brewer colors   # add title
    theme(line = element_blank(),
          axis.text=element_blank(),      # .. tickmarks..
          axis.title=element_blank(),
          legend.position="none",
          plot.margin = unit(c(0, 0, 0, 0), "cm"),
          #legend.text=element_text(size=15),
          #legend.title=element_text(size=17), # .. axis labels..
          panel.background = element_blank(), plot.title = element_text(size=5)) + ggtitle(plot_title)
}
shape_to_ggplot <- function(shape){
  require(broom)
  gg_data <- tidy(shape)
  data <- slot(shape, "data")
  shape[["polyID"]] <- sapply(slot(shape, "polygons"), function(x) slot(x, "ID"))
  gg_data <- merge(gg_data, shape, by.x="id", by.y="polyID")
  return(gg_data)
}
grid_arrange_shared_legend <- function(..., ncol = length(list(...)), nrow = 1, position = c("bottom", "right")) {
  plots <- list(...)
  position <- match.arg(position)
  g <- ggplotGrob(plots[[1]] + theme(legend.position = position))$grobs
  legend <- g[[which(sapply(g, function(x) x$name) == "guide-box")]]
  lheight <- sum(legend$height)
  lwidth <- sum(legend$width)
  gl <- lapply(plots, function(x) x + theme(legend.position="none"))
  gl <- c(gl, ncol = ncol, nrow = nrow)
  combined <- switch(position,
                     "bottom" = arrangeGrob(do.call(arrangeGrob, gl),
                                            legend,
                                            ncol = 1,
                                            heights = unit.c(unit(1, "npc") - lheight, lheight)),
                     "right" = arrangeGrob(do.call(arrangeGrob, gl),
                                           legend,
                                           ncol = 2,
                                           widths = unit.c(unit(1, "npc") - lwidth, lwidth)))
  grid.newpage()
  grid.draw(combined)
  # return gtable invisibly
  invisible(combined)
}
g_legend<-function(a.gplot){
  tmp <- ggplot_gtable(ggplot_build(a.gplot))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  legend <- tmp$grobs[[leg]]
  return(legend)}

Bairro Level Plots

bairro_shp <- readOGR("~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Shapefiles/Shapefiles/Bairro/", "Bairros_from_CTs")

OGR data source with driver: ESRI Shapefile 
Source: "/Users/Sudipta/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Shapefiles/Shapefiles/Bairro", layer: "Bairros_from_CTs"
with 119 features
It has 4 fields

bairro_homs <- read_csv("~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Bairro_SMR_IR_per_yr.csv", 
    col_types = cols(CD_GEOCODB = col_character()))

number of columns of result is not a multiple of vector length (arg 1)17 parsing failures.
row # A tibble: 5 x 5 col     row col   expected              actual file                                                                              expected   <int> <chr> <chr>                 <chr>  <chr>                                                                             actual 1  1854 pop   no trailing characte… e3     '~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Bairro_SMR_I… file 2  1855 pop   no trailing characte… e3     '~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Bairro_SMR_I… row 3  1856 pop   no trailing characte… e3     '~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Bairro_SMR_I… col 4  1857 pop   no trailing characte… e3     '~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Bairro_SMR_I… expected 5  1858 pop   no trailing characte… e3     '~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Bairro_SMR_I…
... ................................. ... ............................................................................................................................ ........ ......................................................................................................................................................................................................................................................... ...... ..................................................................................................................................................... .... ..................................................................................................................................................... ... ..................................................................................................................................................... ... ..................................................................................................................................................... ........ .....................................................................................................................................................
See problems(...) for more details.

bairro_shp_gg <- shape_to_ggplot(bairro_shp)

Regions defined for each Polygons
column name ‘id’ is duplicated in the result

Aggregated plot of counts

bairro_agg_homs <- bairro_homs %>% group_by(CD_GEOCODB) %>% summarize(obs_count=sum(obs_count), SMR=mean(SMR_per_yr), IR=mean(IR_per_yr))

breaks <- quantile(bairro_agg_homs$obs_count, probs=c(seq(0,1,0.15),1))
bairro_shp_gg_count_agg <- bairro_shp_gg[,-1] %>% left_join(bairro_agg_homs[,c(1,2)], by=c("cod_ibge"="CD_GEOCODB")) %>% mutate(cut=cut(obs_count, breaks=breaks, include.lowest=TRUE))

Column `cod_ibge`/`CD_GEOCODB` joining factor and character vector, coercing into character vector

plot_fort_maps2(bairro_shp_gg_count_agg, "cut","Homicide counts", "Aggregate counts for 2001-17 by bairro") + theme(plot.title = element_text(size=10), legend.position = "right")

Aggregated plot of SMRs

breaks <- quantile(bairro_agg_homs$SMR, probs=c(seq(0,1,0.15),1))
bairro_shp_gg_count_agg <- bairro_shp_gg[,-1] %>% left_join(bairro_agg_homs[,c(1,3)], by=c("cod_ibge"="CD_GEOCODB")) %>% mutate(cut=cut(SMR, breaks=breaks, include.lowest=TRUE))

Column `cod_ibge`/`CD_GEOCODB` joining factor and character vector, coercing into character vector

plot_fort_maps2(bairro_shp_gg_count_agg, "cut","Homicide counts", "Aggregate SMR for 2001-17 by bairro") + theme(plot.title = element_text(size=10), legend.position = "right")

Aggregated plot of IRs

breaks <- quantile(bairro_agg_homs$IR, probs=c(seq(0,1,0.15),1))
bairro_shp_gg_count_agg <- bairro_shp_gg[,-1] %>% left_join(bairro_agg_homs[,c(1,4)], by=c("cod_ibge"="CD_GEOCODB")) %>% mutate(cut=cut(IR, breaks=breaks, include.lowest=TRUE))

Column `cod_ibge`/`CD_GEOCODB` joining factor and character vector, coercing into character vector

plot_fort_maps2(bairro_shp_gg_count_agg, "cut","Incidence Rate\n(per 100,000 PYs)", "Aggregate IR for 2001-17 by bairro") + theme(plot.title = element_text(size=10), legend.position = "right")

Count plots by year

breaks <- quantile(bairro_homs$obs_count, probs=c(seq(0,1,0.15),1))
bairro_shp_gg_count <- bairro_shp_gg[,-1] %>% left_join(spread(bairro_homs[,c(1:3)], YOD, obs_count), by=c("cod_ibge"="CD_GEOCODB")) %>% mutate_at(.funs = funs(cut = cut(., breaks=breaks, include.lowest=TRUE)), .vars = vars(`2001`:`2017`))

Column `cod_ibge`/`CD_GEOCODB` joining factor and character vector, coercing into character vector

count_plots <- sapply(2001:2017, function(x) plot_fort_maps2(bairro_shp_gg_count, paste0("`",x,"_cut`"), paste0("Homicide count for ",x), as.character(x)), simplify = FALSE)

library(gridExtra)
library(grid)
library(ggpubr)
library(cowplot)
plot <- count_plots[[1]] + theme(legend.position = "right", legend.text=element_text(size=5), legend.title=element_text(size=5), legend.key.size = unit(0.5,"line")) + guides(fill=guide_legend(title="Homicide Count"))
mylegend <- g_legend(plot)
grid_count <- do.call((plot_grid), c(count_plots, nrow=5)) + draw_grob(mylegend, 0.28, 0.08, 0.28, 0.08) 

IR plots by years

breaks <- quantile(bairro_homs$IR_per_yr, probs=c(seq(0,1,0.15),1), na.rm=TRUE)
bairro_shp_gg_IR <- bairro_shp_gg[,-1] %>% left_join(spread(bairro_homs[,c(1,2,7)], YOD, IR_per_yr), by=c("cod_ibge"="CD_GEOCODB")) %>% mutate_at(.funs = funs(cut = cut(., breaks=breaks, include.lowest=TRUE)), .vars = vars(`2001`:`2017`))

Column `cod_ibge`/`CD_GEOCODB` joining factor and character vector, coercing into character vector

IR_plots <- sapply(2001:2017, function(x) plot_fort_maps2(bairro_shp_gg_IR, paste0("`",x,"_cut`"), paste0("IR for ",x), as.character(x)), simplify = FALSE)
plot <- IR_plots[[1]] + theme(legend.position = "right", legend.text=element_text(size=5), legend.title=element_text(size=5), legend.key.size = unit(0.5,"line")) + guides(fill=guide_legend(title="Homicide IR"))
mylegend <- g_legend(plot)
grid_IR <- do.call((plot_grid), c(count_plots, nrow=5)) + draw_grob(mylegend, 0.28, 0.08, 0.28, 0.08)

SMR plots by years

breaks <- quantile(bairro_homs$SMR_per_yr, probs=c(seq(0,1,0.15),1), na.rm=TRUE)
bairro_shp_gg_SMR <- bairro_shp_gg[,-1] %>% left_join(spread(bairro_homs[,c(1,2,6)], YOD, SMR_per_yr), by=c("cod_ibge"="CD_GEOCODB")) %>% mutate_at(.funs = funs(cut = cut(., breaks=breaks, include.lowest=TRUE)), .vars = vars(`2001`:`2017`))

Column `cod_ibge`/`CD_GEOCODB` joining factor and character vector, coercing into character vector

SMR_plots <- sapply(2001:2017, function(x) plot_fort_maps2(bairro_shp_gg_SMR, paste0("`",x,"_cut`"), paste0("IR for ",x), as.character(x)), simplify = FALSE)
plot <- SMR_plots[[1]] + theme(legend.position = "right", legend.text=element_text(size=5), legend.title=element_text(size=5), legend.key.size = unit(0.5,"line")) + guides(fill=guide_legend(title="Homicide IR"))
mylegend <- g_legend(plot)
grid_SMR <- do.call((plot_grid), c(count_plots, nrow=5)) + draw_grob(mylegend, 0.28, 0.08, 0.28, 0.08)

CT level plots

CT_shp <- readOGR("~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Shapefiles/Shapefiles/CTs/", "Corrected_CTs")

OGR data source with driver: ESRI Shapefile 
Source: "/Users/Sudipta/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/Shapefiles/Shapefiles/CTs", layer: "Corrected_CTs"
with 3044 features
It has 13 fields

CT_homs <- read_csv("~/Documents/Harvard - SM80/Thesis/Fortaleza_Hom_RGit_PRIVATE_Files/CT_SMR_IR_per_yr-agg.csv", 
    col_types = cols(CD_GEOCODI = col_character()))

CT_shp_gg <- shape_to_ggplot(CT_shp)

Regions defined for each Polygons

Some CTs have homicides but non pop. We have to revisit these. For now, I am turning Inf IR and SMR to zero. At most there are 2 homicides per year in these CTs. Also NaN are pop zero, hom zero. Also turning these to zero.

CT_homs <- CT_homs %>% mutate(SMR_per_yr=replace(SMR_per_yr, SMR_per_yr==Inf | is.nan(SMR_per_yr), 0), IR_per_yr=replace(IR_per_yr, IR_per_yr==Inf | is.nan(IR_per_yr), 0)) %>% as.data.frame()

Aggregated plot of counts

CT_agg_homs <- CT_homs %>% group_by(CD_GEOCODI) %>% summarize(obs_count=sum(obs_count), SMR=mean(SMR_per_yr), IR=mean(IR_per_yr))

breaks <- quantile(CT_agg_homs$obs_count, probs=c(seq(0,1,0.15),1))
CT_shp_gg_count_agg <- CT_shp_gg[,-1] %>% left_join(CT_agg_homs[,c(1,2)], by="CD_GEOCODI") %>% mutate(cut=cut(obs_count, breaks=breaks, include.lowest=TRUE))

Column `CD_GEOCODI` joining factor and character vector, coercing into character vector

plot_fort_maps2(CT_shp_gg_count_agg, "cut","Homicide counts", "Aggregate counts for 2001-17 by CT") + theme(plot.title = element_text(size=10), legend.position = "right")

Aggregated plot of SMRs

breaks <- quantile(CT_agg_homs$SMR, probs=c(seq(0,1,0.15),1))
CT_shp_gg_count_agg <- CT_shp_gg[,-1] %>% left_join(CT_agg_homs[,c(1,3)], by="CD_GEOCODI") %>% mutate(cut=cut(SMR, breaks=breaks, include.lowest=TRUE))

Column `CD_GEOCODI` joining factor and character vector, coercing into character vector

plot_fort_maps2(CT_shp_gg_count_agg, "cut","SMR", "Aggregate SMR for 2001-17 by CT") + theme(plot.title = element_text(size=10), legend.position = "right")

Aggregated plot of IRs

breaks <- quantile(CT_agg_homs$IR, probs=c(seq(0,1,0.15),1))
CT_shp_gg_count_agg <- CT_shp_gg[,-1] %>% left_join(CT_agg_homs[,c(1,4)], by="CD_GEOCODI") %>% mutate(cut=cut(IR, breaks=breaks, include.lowest=TRUE))

Column `CD_GEOCODI` joining factor and character vector, coercing into character vector

plot_fort_maps2(CT_shp_gg_count_agg, "cut","Incidence Rate\n(per 100,000 PYs)", "Aggregate IR for 2001-17 by CT") + theme(plot.title = element_text(size=10), legend.position = "right")

Count plots by year

breaks <- c(0,1,2,3,4,10,35)
CT_shp_gg_count <- CT_shp_gg[,-1] %>% left_join(spread(CT_homs[,c(1:3)], YOD, obs_count), by="CD_GEOCODI") %>% mutate_at(.funs = funs(cut = cut(., breaks=breaks, include.lowest=TRUE)), .vars = vars(`2001`:`2017`))

Column `CD_GEOCODI` joining factor and character vector, coercing into character vector

count_plots <- sapply(2001:2017, function(x) plot_fort_maps2(CT_shp_gg_count, paste0("`",x,"_cut`"), paste0("Homicide count for ",x), as.character(x), size=0.05), simplify = FALSE)

library(gridExtra)
library(grid)
library(ggpubr)
library(cowplot)
plot <- count_plots[[1]] + theme(legend.position = "right", legend.text=element_text(size=5), legend.title=element_text(size=5), legend.key.size = unit(0.5,"line")) + guides(fill=guide_legend(title="Homicide Count"))
mylegend <- g_legend(plot)
grid_count <- do.call((plot_grid), c(count_plots, nrow=5)) + draw_grob(mylegend, 0.28, 0.08, 0.28, 0.08) 

grid_count

IR plots by years

breaks <- c(0,5,10,15,50,100,200,15000)
CT_shp_gg_IR <- CT_shp_gg[,-1] %>% left_join(spread(CT_homs[,c(1,2,7)], YOD, IR_per_yr), by="CD_GEOCODI") %>% mutate_at(.funs = funs(cut = cut(., breaks=breaks, labels=(c("[0,5]","(5,10]","(10,15]","(15,50]","(50,100]","(100,200]",">200")), include.lowest=TRUE)), .vars = vars(`2001`:`2017`))

Column `CD_GEOCODI` joining factor and character vector, coercing into character vector

IR_plots <- sapply(2001:2017, function(x) plot_fort_maps2(CT_shp_gg_IR, paste0("`",x,"_cut`"), paste0("IR for ",x), as.character(x)), simplify = FALSE)
plot <- IR_plots[[1]] + theme(legend.position = "right", legend.text=element_text(size=5), legend.title=element_text(size=5), legend.key.size = unit(0.5,"line")) + guides(fill=guide_legend(title="Homicide IR"))
mylegend <- g_legend(plot)
grid_IR <- do.call((plot_grid), c(count_plots, nrow=5)) + draw_grob(mylegend, 0.28, 0.08, 0.28, 0.08)

grid_IR

SMR plots by years

breaks <- c(0,1,2,5,10,15,50,400)
CT_shp_gg_SMR <- CT_shp_gg[,-1] %>% left_join(spread(CT_homs[,c(1,2,6)], YOD, SMR_per_yr), by="CD_GEOCODI") %>% mutate_at(.funs = funs(cut = cut(., breaks=breaks, labels=(c("[0,1]","(1,2]","(2,5]","(5,10]","(10,15]","(15,50]",">50")), include.lowest=TRUE)), .vars = vars(`2001`:`2017`))

Column `CD_GEOCODI` joining factor and character vector, coercing into character vector

SMR_plots <- sapply(2001:2017, function(x) plot_fort_maps2(CT_shp_gg_SMR, paste0("`",x,"_cut`"), paste0("IR for ",x), as.character(x)), simplify = FALSE)
plot <- SMR_plots[[1]] + theme(legend.position = "right", legend.text=element_text(size=5), legend.title=element_text(size=5), legend.key.size = unit(0.5,"line")) + guides(fill=guide_legend(title="SMR"))
mylegend <- g_legend(plot)
grid_SMR <- do.call((plot_grid), c(count_plots, nrow=5)) + draw_grob(mylegend, 0.28, 0.08, 0.28, 0.08)

grid_SMR

CT_shp_gg_SMR$SMR_time1 <- CT_shp_gg_SMR$`2001` + CT_shp_gg_SMR$`2002` + CT_shp_gg_SMR$`2003` + CT_shp_gg_SMR$`2004` + CT_shp_gg_SMR$`2005` + CT_shp_gg_SMR$`2006`
CT_shp_gg_SMR$SMR_time2 <- CT_shp_gg_SMR$`2001` + CT_shp_gg_SMR$`2002` + CT_shp_gg_SMR$`2003` + CT_shp_gg_SMR$`2004` + CT_shp_gg_SMR$`2005` + CT_shp_gg_SMR$`2006`